fig3
For this analysis, we used the following parameters:
pars <- as.data.frame(do.call("rbind", params))
colnames(pars) <- "parameters_used"
parsLoad Seurat object.
git_path <- params$git_path
figure_out_path <- params$figure_out_path
sobj <- readRDS(file=file.path(git_path,params$object))Load DESeq2 pseudobulk results. For full documentation of DESeq2 analysis, refer to Scritps/DESeq2.Rmd.
res_df <- read.delim(file.path(git_path,"Figures/data/de_results_pseudoid_analysis.tsv"))
selected_contrasts <- c("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx",
"26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx",
"26w_SPF_5days_fx_vs_26w_SPF_2days_fx",
"26w_nonSPF_2days_fx_vs_12w_nonSPF_2days_fx",
"26w_SPF_2days_fx_vs_12w_nonSPF_2days_fx",
"26w_SPF_2days_fx_vs_26w_nonSPF_2days_fx",
"12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx",
"12w_nonSPF_5days_proximal_vs_12w_nonSPF_5days_fx",
"12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_proximal",
"26w_nonSPF_5days_distal_vs_26w_nonSPF_5days_fx",
"26w_nonSPF_5days_proximal_vs_26w_nonSPF_5days_fx",
"26w_nonSPF_5days_distal_vs_26w_nonSPF_5days_proximal",
"26w_SPF_5days_distal_vs_26w_SPF_5days_fx",
"26w_SPF_5days_proximal_vs_26w_SPF_5days_fx",
"26w_SPF_5days_distal_vs_26w_SPF_5days_proximal")
cs <- data.frame(contrast=unique(res_df$contrast),
new=selected_contrasts)
res_df <- res_df %>%
left_join(cs) %>%
dplyr::rename(old_contrast=contrast,
contrast=new)
gene_info <- read.delim(file.path(git_path,params$gene_info), header=F)
colnames(gene_info) <- c("gene_id","gene_name","gene_type")
res_df <- res_df %>%
left_join(gene_info)for (comparison in unique(res_df$contrast)){
p <- res_df %>%
dplyr::filter(!is.na(padj), contrast == comparison) %>%
ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .05)) +
geom_point(size=.5, shape=1) +
geom_label_repel(data=res_df %>%
dplyr::filter(contrast == comparison,padj < .05),
aes(label=gene_name),
segment.size=.25,segment.alpha=.5,size=3,color='black',max.overlaps=40,
label.padding=0.05,label.size=NA) +
scale_color_manual(values=c('black','red')) +
theme_classic() +
theme(legend.position='none',
strip.background=element_blank())+
ggtitle(comparison)+
facet_wrap(~cell_type, ncol = 2, scales = "free")
print(p)
}
res_df %>%
dplyr::arrange(padj) %>%
dplyr::filter(!is.na(padj) & (padj < .05)) %>%
dplyr::mutate_if(is.numeric, signif, 2) %>%
DT::datatable(rownames=FALSE, extensions = 'Buttons', options = list(dom='frtBip', buttons = c('csv')))library(tmod)
# this is to document how the mouse orthologues were mapped for the tmod pathways
# data(tmod)
#
# library(orthomapper)
# mousehum <- orthologs(10090, 9606)
# mousehum$h_symbol <- entrez_annotate(mousehum$T.9606)$SYMBOL
# mousehum$m_symbol <- entrez_annotate(mousehum$T.10090)$SYMBOL
#
# mousehum <- mousehum[order(mousehum$h_symbol),]
# rows_tmod <- which(mousehum$h_symbol %in% tmod$gv)
# mgenes <- mousehum$m_symbol[rows_tmod]
#
# list_tmod <- sapply(tmod$gs2gv, function(x) tmod$gv[x] )
# list_tmod <- sapply(list_tmod, function(x) x[order(x)])
# list_tmod <- sapply(list_tmod, function(x) x[which(x %in% mousehum$h_symbol)])
# list_tmod <- sapply(list_tmod, function(x) mousehum[which(mousehum$h_symbol %in% x), "m_symbol"])
# list_tmod <- sapply(list_tmod, function(x) which(mgenes %in% x))
#
# tmod_mm <- tmod
#
# tmod_mm$gv <- mgenes
# tmod_mm$gs2gv <- list_tmod
# save tmod mouse object
# saveRDS(tmod_mm, "tmod_mm.rds")
tmod_mm <- readRDS(file.path(git_path,"Figures/data/tmod_mm.rds"))
df2tmod <- function(df, gene_id_col=ncol(df), module_id_col=1, module_title_col=2) {
require(tmod, quietly=TRUE)
gene_ids <- df[[gene_id_col]]
module_ids <- df[[module_id_col]]
m2g <- tapply(gene_ids, module_ids, list)
df[[gene_id_col]] <- NULL
df <- df[!duplicated(df[[module_id_col]]), ]
colnames(df)[module_id_col] <- "ID"
colnames(df)[module_title_col] <- "Title"
makeTmod(modules=df, modules2genes=m2g)
}
df <- as.data.frame(msigdbr::msigdbr(species='Mus musculus'))
df <- df[ , c("gs_name", "gs_id", "gs_cat", "gs_subcat", "gene_symbol") ]
colnames(df) <- c("Title", "ID", "Category", "Subcategory", "GeneID")
msig <- df2tmod(df[!is.na(df$GeneID),], gene_id_col=ncol(df), module_id_col=2, module_title_col=1)
sel <- (msig$gs$Category %in% c("H")) | (msig$gs$Subcategory %in% c('CP:REACTOME','GO:BP','CP:KEGG'))res <- res_df %>%
filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)
genes <- lapply(res, function(x) x %>%
filter(!is.na(padj)) %>%
arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]
lgenes <- lapply(genes_reduced, function(x) x$gene_name)
res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)
res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))Figure 3A
#################### Volcano plot #######################
DF <- res_df[grepl("26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP"
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange < (-0.5) & DF$padj < 0.1)] <- "down"
(p12 <- DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(!is.na(padj)) %>%
ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) +
geom_point(size=.5, shape=1) +
geom_label_repel(data=DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(direction=="up"|direction=="down") ,
aes(label=gene_name),
segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
label.padding=0.05,label.size=NA) +
scale_color_manual(values=c("blue","grey","red"))+ theme_classic() +
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8,hjust = 0.5,),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
axis.ticks.x=element_blank(),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_blank()) +
labs(title = "",x = "log2FoldChange (5 days vs. 2 days)") +
ggtitle("Immune-aged: hematoma"))
Figure S3A
#################### Volcano plot #######################
DF <- res_df[grepl("26w_SPF_5days_fx_vs_26w_SPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP"
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange < (-0.5) & DF$padj < 0.1)] <- "down"
(sf3a <- DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(!is.na(padj)) %>%
ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) +
geom_point(size=.5, shape=1) +
geom_label_repel(data=DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(direction=="up"|direction=="down") ,
aes(label=gene_name),
segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
label.padding=0.05,label.size=NA) +
scale_color_manual(values=c("blue","grey","red"))+ theme_classic() +
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8,hjust = 0.5,),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
axis.ticks.x=element_blank(),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_blank()) +
labs(title = "",x = "log2FoldChange (5 days vs. 2 days)")+
ggtitle("Aged: hematoma"))
Figure S3B
#################### Volcano plot #######################
DF <- res_df[grepl("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", res_df$contrast), ]
# add a column of NAs
DF$direction <- "N"
# if log2Foldchange > 0.5 and pvalue < 0.05, set as "UP"
DF$direction[(DF$log2FoldChange > 0.5 & DF$padj < 0.1)] <- "up"
DF$direction[(DF$log2FoldChange < (-0.5) & DF$padj < 0.1)] <- "down"
(sf3b <- DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(!is.na(padj)) %>%
ggplot(aes(x=log2FoldChange,y=-log10(pvalue),col=direction)) +
geom_point(size=.5, shape=1) +
geom_label_repel(data=DF %>%
dplyr::filter(gene_type=="protein_coding") %>%
dplyr::filter(direction=="up"|direction=="down") ,
aes(label=gene_name),
segment.size=.25,segment.alpha=.5,size=2.5,color='black',max.overlaps=30,
label.padding=0.05,label.size=NA) +
scale_color_manual(values=c("blue","grey","red"))+ theme_classic() +
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8,hjust = 0.5,),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
axis.ticks.x=element_blank(),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_blank()) +
labs(title = "",x = "log2FoldChange (5 days vs. 2 days)") +
ggtitle("Young: hematoma"))
Panel plot for the immune-aged hematoma.
res <- res_df %>%
filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)
genes <- lapply(res, function(x) x %>%
filter(!is.na(padj)) %>%
arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]
lgenes <- lapply(genes_reduced, function(x) x$gene_name)
res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)
res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))
x <- genes_reduced$"26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)
names(sgenes) <- "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"
p10 <- ggPanelplot(res.tmod.filtered, sgenes = sgenes,q_cutoff = 1e-7,cluster=F) +
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8,hjust = 0.5,),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
axis.ticks.x=element_blank(),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_blank()) +
labs(title = "Immune-aged: \nhematoma, \n5 days vs 2 days")
p10
For young and aged no enrichment is found, so no figures are generated here.
res <- res_df %>%
filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)
genes <- lapply(res, function(x) x %>%
filter(!is.na(padj)) %>%
arrange(pvalue))
genes_reduced <- genes["12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"]
lgenes <- lapply(genes_reduced, function(x) x$gene_name)
res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)
res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))
x <- genes_reduced$"12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)
names(sgenes) <- "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
# sf3c <- ggPanelplot(res.tmod.filtered, sgenes = sgenes, q_cutoff = 1e-7,cluster=F) +
# theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8,hjust = 0.5,),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# axis.ticks.x=element_blank(),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# strip.text.x = element_blank()) +
# labs(title = "Young: \nhematoma, \n5 days vs 2 days")
# sf3cres <- res_df %>%
filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)
genes <- lapply(res, function(x) x %>%
filter(!is.na(padj)) %>%
arrange(pvalue))
genes_reduced <- genes["26w_SPF_5days_fx_vs_26w_SPF_2days_fx"]
lgenes <- lapply(genes_reduced, function(x) x$gene_name)
res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)
res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))
x <- genes_reduced$"26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)
names(sgenes) <- "12w_nonSPF_5days_distal_vs_12w_nonSPF_5days_fx"
# sf3d <- ggPanelplot(res.tmod.filtered, sgenes = sgenes, q_cutoff = 1e-7,cluster=F) +
# theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8,hjust = 0.5,),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# axis.ticks.x=element_blank(),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# strip.text.x = element_blank()) +
# labs(title = "Aged: \nhematoma, \n5 days vs 2 days")
# sf3dGenerating figure 4C
res <- res_df %>%
filter(gene_type=="protein_coding")
res <- split(res, f = res$contrast)
genes <- lapply(res, function(x) x %>%
filter(!is.na(padj)) %>%
arrange(pvalue))
genes_reduced <- genes["26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"]
lgenes <- lapply(genes_reduced, function(x) x$gene_name)
res.tmod <- lapply(lgenes, tmodCERNOtest, mset = tmod_mm)
res.tmod.filtered <- lapply(res.tmod, function(x) subset(x, adj.P.Val < 1e-3 & AUC > .75))
v <- unlist(lgenes, use.names=FALSE)#2 for non SPF
down<-x %>% filter(log2FoldChange < 0.5, padj<0.1) %>% pull(gene_name)
modmetabo<-tmod_mm
modmetabo$gs[modmetabo$gs$ID %in% c("LI.M47.0"), ]path_positions<- modmetabo$gs2gv[[which(modmetabo$gs$ID=="LI.M47.0")]]
genes<-modmetabo$gv[path_positions]
Genes_downs = down[down %in% genes]
Genes_all = v[v %in% genes]
gene.colors <- c("grey","grey","blue" ,"grey","blue" ,"grey", "blue", "grey","grey","grey","grey", "grey","grey", "grey","grey","grey","grey", "grey", "grey", "grey","grey","grey","grey","grey","grey","grey" ,"grey","grey","grey","blue","grey" ,"grey","grey","grey","grey")
names(gene.colors)<- genes
######### Evidence plot ###############################
update_geom_defaults("text_repel", list(size = 2.5))
(p11 <- ggEvidencePlot(v, "LI.M47.0", mset = tmod_mm, gene.colors = gene.colors,filter=F) +
theme_classic() +
theme(text = element_text(family = "Arial",size=8),
plot.title = element_text(size = 8,hjust = 0.5,),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
axis.ticks.x=element_blank(),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
legend.position = "none",
strip.text.x = element_blank()) + labs(title = "Genes in LI.M47.0"))
Generating
Figure E & F
############# Downsample ################################
# Setting seeds
set.seed(10)
Idents(sobj) <- "group"
Downsample<- subset(sobj, downsample=500)
############# Umap plot ################################
modmetabo<-tmod_mm
modmetabo$gs[modmetabo$gs$ID %in% c("LI.M47.0"), ]path_positions<- modmetabo$gs2gv[[which(modmetabo$gs$ID=="LI.M47.0")]]
genes<-modmetabo$gv[path_positions]
Genes_Bcells= v[v %in% genes]
DefaultAssay(Downsample) <- "RNA"
Downsample <- AddModuleScore(Downsample,
features = list(Genes_Bcells),
name="Genes_Bcells")
Idents(Downsample)<-Downsample$timepoint
testdata <-subset(Downsample, idents = c("2days"))
Idents(testdata)<-testdata$compartments
test <-subset(testdata, idents = c("Hematoma"))
Idents(test)<-test$Age_raw
Immunoaged <- subset(test, idents = c("Immunoaged"))
Idents(Immunoaged)<-Immunoaged$CellTypes
# Plot module scores bcells 2 days
p13<-FeaturePlot(Immunoaged,
features = "Genes_Bcells1",label = FALSE, repel = TRUE) +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
guides(colour=guide_colorbar(title="Gene module \nscore"))+
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_text(size = 8),
legend.position = "top",
legend.justification = "center") + labs(title = "Immune-aged: hematoma, 2 days") & NoAxes()
p13
################################# 5 days ############################################
Idents(Downsample)<-Downsample$timepoint
testdata <-subset(Downsample, idents = c("5days"))
Idents(testdata)<-testdata$compartments
test <-subset(testdata, idents = c("Hematoma"))
Idents(test)<-test$Age_raw
Immunoaged <- subset(test, idents = c("Immunoaged"))
Idents(Immunoaged)<-Immunoaged$CellTypes
# Plot module scores bcells 5 days
p14<-FeaturePlot(Immunoaged,
features = "Genes_Bcells1",label = FALSE, repel = TRUE) +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
guides(colour=guide_colorbar(title="Gene module \nscore"))+
theme(text = element_text(family = "Arial"),
plot.title = element_text(size = 8),
axis.title = element_text(size = 8),
axis.text = element_text(size = 8),
legend.text = element_text(size =8),
legend.title = element_text(size = 8),
strip.text.x = element_text(size = 8),
legend.position = "top",
legend.justification = "center") + labs(title = "Immune-aged: hematoma, 5 days") & NoAxes()
p14
Interferon - which cell type? UMAP
# ############# Downsample ################################
# # Setting seeds
# set.seed(10)
# Idents(sobj) <- "group"
# Downsample<- subset(sobj, downsample=500)
# ############# Umap plot ################################
#
# ifn <- read.delim(file.path(git_path,"Figures/data/signatures_mm_signatuR.txt"), header=T) %>%
# filter(signature=="IfnResp") %>%
# pull(gene)
#
# genes <- ifn
#
# Genes_Bcells= v[v %in% genes]
#
# DefaultAssay(Downsample) <- "RNA"
# Downsample <- AddModuleScore(Downsample,
# features = list(Genes_Bcells),
# name="IfnResponse")
#
# Idents(Downsample)<-Downsample$timepoint
# testdata <-subset(Downsample, idents = c("2days"))
# Idents(testdata)<-testdata$compartments
# test <-subset(testdata, idents = c("Hematoma"))
# Idents(test)<-test$Age_raw
# Immunoaged <- subset(test, idents = c("Immunoaged"))
# Idents(Immunoaged)<-Immunoaged$CellTypes
# # Plot module scores bcells 2 days
# p13<-FeaturePlot(Immunoaged,
# features = "IfnResponse1",label = FALSE, repel = TRUE) +
# scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
# guides(colour=guide_colorbar(title="Gene module \nscore"))+
# theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# strip.text.x = element_text(size = 8),
# legend.position = "top",
# legend.justification = "center") + labs(title = "Immune-aged: hematoma, 2 days") & NoAxes()
#
# p13
#
# ################################# 5 days ############################################
# Idents(Downsample)<-Downsample$timepoint
# testdata <-subset(Downsample, idents = c("5days"))
# Idents(testdata)<-testdata$compartments
# test <-subset(testdata, idents = c("Hematoma"))
# Idents(test)<-test$Age_raw
# Immunoaged <- subset(test, idents = c("Immunoaged"))
# Idents(Immunoaged)<-Immunoaged$CellTypes
# # Plot module scores bcells 5 days
# p14<-FeaturePlot(Immunoaged,
# features = "IfnResponse1",label = FALSE, repel = TRUE) +
# scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")), limits = c(0.0, 1.5)) +
# guides(colour=guide_colorbar(title="Gene module \nscore"))+
# theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# strip.text.x = element_text(size = 8),
# legend.position = "top",
# legend.justification = "center") + labs(title = "Immune-aged: hematoma, 5 days") & NoAxes()
#
# p14Generating Figure 4 D
library(ggpubr)
my_comparisons <- list(c("2days", "5days"))
DefaultAssay(Downsample) <- "RNA"
p15<-Downsample[[]] %>%
subset(compartments == "Hematoma") %>%
subset(Age_raw=="Immunoaged") %>%
filter(CellTypes %in% c("B Cells")) %>%
select(sample,group, Age_raw,Genes_Bcells1,compartments,timepoint,CellTypes) %>%
#mutate(Age_raw=factor(Age_raw, levels = c("Young","Aged","ImmunoAged"))) %>%
ggplot(aes(timepoint, Genes_Bcells1))+
geom_boxplot(size = .4, outlier.size = 0) +
geom_point(aes(color = timepoint), size = .6) +
#facet_wrap(~Age_raw, scales="free_x")+
ylim(0, 2) +
stat_compare_means(comparisons = my_comparisons,label="p.signif",
method = "t.test",label.y = 1.5, label.x = "5days",size = 6.1167,hide.ns = TRUE) + theme_classic() +
theme(axis.title = element_text(size = 8),
axis.text = element_text(size = 8, family = "Arial"),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_text(size =8),
legend.text = element_text(size =8),
legend.title = element_text(size =8),
strip.text.x = element_text(size = 8),
strip.text.y = element_text(size = 8),
plot.title = element_text(size =8,hjust = 0.5,))+
scale_color_manual(name= "Time",values = c("dodgerblue3","orange2")) + labs(title ="Immune-aged: hematoma", x = c(color =''), y = 'Gene module score of LI.M47.0')
p15
Generating Figure 3
(wrap_elements(p12) + wrap_elements(p10)) /
(wrap_elements(p11) + wrap_elements(p15) ) /
(wrap_elements(p13) + wrap_elements(p14) ) +
plot_layout(heights=c(1,.8,1.1)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 8, face = "bold",hjust = -.5, vjust = 0))
if (file.exists(file.path(figure_out_path,"fig3.pdf"))){
try(ggsave(file.path(figure_out_path,"fig3.pdf"), device = cairo_pdf,width=8,height=9), silent=T)
warning("file not written, please close file in the other application and rerun.")
} else {
ggsave(file.path(figure_out_path,"fig3.pdf"), device = cairo_pdf,width=8,height=9)
}Generating Suppl Figure 3
(wrap_elements(sf3a) + wrap_elements(sf3b)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 8, face = "bold",hjust = -.5, vjust = 0))
if (file.exists(file.path(figure_out_path,"SFig3.pdf"))){
try(ggsave(file.path(figure_out_path,"SFig3.pdf"), device = cairo_pdf,width=8,height=9), silent=T)
warning("file not written, please close file in the other application and rerun.")
} else {
ggsave(file.path(figure_out_path,"SFig3.pdf"), device = cairo_pdf,width=8,height=3)
}supplement figure (SFig A to D)
# library(RColorBrewer)
# #################### Volcano plot #######################
# DF <- res_df[grepl("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", res_df$contrast), ]
# (SFig4_1 <- DF %>%
# dplyr::filter(!is.na(padj)) %>%
# ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .1)) +
# geom_point(size=.5, shape=1) +
# geom_label_repel(data=DF %>%
# dplyr::filter(padj < .1),
# aes(label=gene_name),
# segment.size=.25,segment.alpha=.5,size=4,color='black',max.overlaps=15,
# label.padding=0.05,label.size=NA) +
# scale_color_manual(values=c('black','red'))+ theme_classic() + theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8,hjust = 0.5,),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# axis.ticks.x=element_blank(),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# legend.position = "none",
# strip.text.x = element_blank()) + labs(title = "Young",x = "log2FoldChange (5 days vs. 2 days)"))
#
# DF <- res_df[grepl("26w_SPF_5days_fx_vs_26w_SPF_2days_fx", res_df$contrast), ]
# (SFig4_2 <- DF %>%
# dplyr::filter(!is.na(padj)) %>%
# ggplot(aes(x=log2FoldChange,y=-log10(pvalue),color=padj < .1)) +
# geom_point(size=.5, shape=1) +
# geom_label_repel(data=DF %>%
# dplyr::filter(padj < .1),
# aes(label=gene_name),
# segment.size=.25,segment.alpha=.5,size=4,color='black',max.overlaps=15,
# label.padding=0.05,label.size=NA) +
# scale_color_manual(values=c('black','red'))+ theme_classic() + theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8,hjust = 0.5,),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# axis.ticks.x=element_blank(),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# legend.position = "none",
# strip.text.x = element_blank()) + labs(title = "Aged",x = "log2FoldChange (5 days vs. 2 days)"))Generating supplement Figure E
# library(tmod)
# library(disco)
# dataset_1 <- res_df[res_df$contrast == "12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", ]
# dataset_2 <- res_df[res_df$contrast == "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", ]
# dataset_3 <- res_df[res_df$contrast == "26w_SPF_5days_fx_vs_26w_SPF_2days_fx", ]
# names<- c("Aged","Immunoaged")
# orthologs <- NULL
# orthologs <- makeMatchedOrtholog(names,dataset_3$gene_name, dataset_3$log2FoldChange,dataset_2$log2FoldChange, dataset_3$padj,dataset_2$padj,row.names = NULL, extra = NULL)
# #ds <- disco.score(orthologs)
# #plotDisco(orthologs, ds)
#
# both<- orthologs$lfc %>% as.data.frame()
# both$gene_name <- orthologs$genes
# genes <- res_df[res_df$contrast == c("26w_SPF_5days_fx_vs_26w_SPF_2days_fx","26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"), ]
# Genes <- genes[genes$padj < 0.05 &
# !is.na(genes$padj) &
# abs(genes$log2FoldChange) > 0.5,] %>% pull(gene_name)Plot logfoldchange against one another
# # add a column of NAs
# both$colours <- "NO"
# # if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP"
# both$colours[(both[,1] > 0.5 & dataset_3$padj< 0.05)] <- "UP"
# both$colours[(both[,2] > 0.5 & dataset_2$padj< 0.05)] <- "UP"
# # if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
# both$colours[(both[,1] < -0.5 & dataset_3$padj< 0.05)] <- "DOWN"
# both$colours[(both[,2] < -0.5 & dataset_2$padj< 0.05)] <- "DOWN"
# both$gene_set_labels <- ifelse(both$gene_name %in% Genes,both$gene_name, NA)# SFig4_5 <- ggplot(both, aes(both[,1], both[,2],colour = colours)) + geom_point(shape=1, size = .5) + geom_text_repel(aes(label = gene_set_labels), size = 3, max.overlaps = 100) + ggtitle("Hematoma: Aged vs Immunoaged") + xlab("log2Foldchange (Aged)") + ylab("log2Foldchange (Immunoaged)") + scale_color_manual(values=c("blue","grey","red")) + theme_classic() + theme(plot.title = element_text(size = 8, hjust=0.5))
# SFig4_5Supplemenyt Figure F: disco plot young vs immunoaged
# library(tmod)
# library(disco)
# dataset_1 <- res_df[res_df$contrast == "12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx", ]
# dataset_2 <- res_df[res_df$contrast == "26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx", ]
# dataset_3 <- res_df[res_df$contrast == "26w_SPF_5days_fx_vs_26w_SPF_2days_fx", ]
# names<- c("Young","Immunoaged")
# orthologs <- NULL
# orthologs <- makeMatchedOrtholog(names,dataset_1$gene_name, dataset_1$log2FoldChange,dataset_2$log2FoldChange, dataset_1$padj,dataset_2$padj,row.names = NULL, extra = NULL)
# #ds <- disco.score(orthologs)
# #plotDisco(orthologs, ds)
#
# both<- orthologs$lfc %>% as.data.frame()
# both$gene_name <- orthologs$genes
# genes <- res_df[res_df$contrast == c("12w_nonSPF_5days_fx_vs_12w_nonSPF_2days_fx","26w_nonSPF_5days_fx_vs_26w_nonSPF_2days_fx"), ]
# Genes <- genes[genes$padj < 0.05 &
# !is.na(genes$padj) &
# abs(genes$log2FoldChange) > 0.5,] %>% pull(gene_name)Plot logfoldchange against one another
# # add a column of NAs
# both$colours <- "NO"
# # if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP"
# both$colours[(both[,1] > 0.5 & dataset_1$padj< 0.05)] <- "UP"
# both$colours[(both[,2] > 0.5 & dataset_2$padj< 0.05)] <- "UP"
# # if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
# both$colours[(both[,1] < -0.5 & dataset_1$padj< 0.05)] <- "DOWN"
# both$colours[(both[,2] < -0.5 & dataset_2$padj< 0.05)] <- "DOWN"
# both$gene_set_labels <- ifelse(both$gene_name %in% Genes,both$gene_name, NA)# SFig4_6 <- ggplot(both, aes(both[,1], both[,2],colour = colours)) + geom_point(shape=1, size = .5) + geom_text_repel(aes(label = gene_set_labels), size = 3, max.overlaps = 100) + ggtitle("Hematoma:Young vs Immunoaged") + xlab("log2Foldchange (Young)") + ylab("log2Foldchange (Immunoaged)") + scale_color_manual(values=c("blue","grey","red")) + theme_classic() + theme(plot.title = element_text(size = 8, hjust=0.5))
# SFig4_6generating tmod plot for figure 4
# x <- genes$"26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
# sgenes <- tmodDecideTests(g = x$gene_name, lfc = x$log2FoldChange, pval = x$padj, lfc.thr = 0.5, pval.thr = 0.1, mset = tmod_mm)
#
# names(sgenes) <- "26w_SPF_5days_fx_vs_26w_SPF_2days_fx"
# p10 <- ggPanelplot(res.tmod.filtered, sgenes = sgenes,q_cutoff = 1e-12) +
# theme(text = element_text(family = "Arial"),
# plot.title = element_text(size = 8,hjust = 0.5,),
# axis.title = element_text(size = 8),
# axis.text = element_text(size = 8),
# axis.ticks.x=element_blank(),
# legend.text = element_text(size =8),
# legend.title = element_text(size = 8),
# strip.text.x = element_blank()) +
# labs(title = "Aged: \nhematoma, \n5 days vs 2 days")
# p10Generating supplement Figure 4
# library(patchwork)
# SFig4_1+ SFig4_2 + SFig4_6 + SFig4_5 +
# plot_annotation(tag_levels = 'A') + plot_layout(ncol = 2)
# ggsave(file.path(figure_out_path,"SFig4.pdf"), device = cairo_pdf, width=8,height=9)sessionInfo()## R version 4.3.3 (2024-02-29)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/milekm/work/miniconda3/envs/monocle/lib/libopenblasp-r0.3.27.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Berlin
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggpubr_0.6.0 tmod_0.50.13
## [3] patchwork_1.2.0 RColorBrewer_1.1-3
## [5] stringr_1.5.1 ggrepel_0.9.5
## [7] DESeq2_1.42.0 SummarizedExperiment_1.32.0
## [9] Biobase_2.62.0 MatrixGenerics_1.14.0
## [11] matrixStats_1.3.0 GenomicRanges_1.54.1
## [13] GenomeInfoDb_1.38.1 IRanges_2.36.0
## [15] S4Vectors_0.40.2 BiocGenerics_0.48.1
## [17] cowplot_1.1.3 ggplot2_3.5.1
## [19] gtools_3.9.5 tidyr_1.3.1
## [21] dplyr_1.1.4 Seurat_5.1.0
## [23] SeuratObject_5.0.2 sp_2.1-4
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.22 splines_4.3.3 later_1.3.2
## [4] plotwidgets_0.5.1 bitops_1.0-7 tibble_3.2.1
## [7] polyclip_1.10-7 XML_3.99-0.17 fastDummies_1.7.3
## [10] lifecycle_1.0.4 rstatix_0.7.2 globals_0.16.3
## [13] lattice_0.22-6 MASS_7.3-60.0.1 crosstalk_1.2.1
## [16] backports_1.5.0 magrittr_2.0.3 plotly_4.10.4
## [19] sass_0.4.9 rmarkdown_2.27 jquerylib_0.1.4
## [22] yaml_2.3.9 httpuv_1.6.15 sctransform_0.4.1
## [25] spam_2.10-0 spatstat.sparse_3.1-0 reticulate_1.38.0
## [28] pbapply_1.7-2 abind_1.4-5 zlibbioc_1.48.0
## [31] Rtsne_0.17 purrr_1.0.2 msigdbr_7.5.1
## [34] RCurl_1.98-1.16 GenomeInfoDbData_1.2.11 irlba_2.3.5.1
## [37] listenv_0.9.1 spatstat.utils_3.0-5 pheatmap_1.0.12
## [40] goftest_1.2-3 RSpectra_0.16-2 spatstat.random_3.3-1
## [43] fitdistrplus_1.2-1 parallelly_1.37.1 leiden_0.4.3.1
## [46] codetools_0.2-20 DelayedArray_0.28.0 DT_0.33
## [49] tidyselect_1.2.1 farver_2.1.2 spatstat.explore_3.3-1
## [52] jsonlite_1.8.8 progressr_0.14.0 ggridges_0.5.6
## [55] survival_3.7-0 tools_4.3.3 tagcloud_0.6
## [58] ica_1.0-3 Rcpp_1.0.13 glue_1.7.0
## [61] gridExtra_2.3 SparseArray_1.2.2 xfun_0.46
## [64] withr_3.0.0 fastmap_1.2.0 fansi_1.0.6
## [67] caTools_1.18.2 digest_0.6.36 R6_2.5.1
## [70] mime_0.12 colorspace_2.1-0 scattermore_1.2
## [73] tensor_1.5 spatstat.data_3.1-2 utf8_1.2.4
## [76] generics_0.1.3 data.table_1.15.4 httr_1.4.7
## [79] htmlwidgets_1.6.4 S4Arrays_1.2.0 uwot_0.2.2
## [82] pkgconfig_2.0.3 gtable_0.3.5 lmtest_0.9-40
## [85] XVector_0.42.0 htmltools_0.5.8.1 carData_3.0-5
## [88] dotCall64_1.1-1 scales_1.3.0 png_0.1-8
## [91] spatstat.univar_3.0-0 knitr_1.48 reshape2_1.4.4
## [94] nlme_3.1-165 cachem_1.1.0 zoo_1.8-12
## [97] KernSmooth_2.23-24 parallel_4.3.3 miniUI_0.1.1.1
## [100] pillar_1.9.0 grid_4.3.3 vctrs_0.6.5
## [103] RANN_2.6.1 gplots_3.1.3.1 promises_1.3.0
## [106] car_3.1-2 xtable_1.8-4 cluster_2.1.6
## [109] beeswarm_0.4.0 evaluate_0.24.0 cli_3.6.3
## [112] locfit_1.5-9.9 compiler_4.3.3 rlang_1.1.4
## [115] crayon_1.5.3 ggsignif_0.6.4 future.apply_1.11.2
## [118] labeling_0.4.3 plyr_1.8.9 stringi_1.8.4
## [121] viridisLite_0.4.2 deldir_2.0-4 BiocParallel_1.36.0
## [124] babelgene_22.9 munsell_0.5.1 lazyeval_0.2.2
## [127] spatstat.geom_3.3-2 Matrix_1.6-5 RcppHNSW_0.6.0
## [130] future_1.33.2 shiny_1.8.1.1 highr_0.11
## [133] ROCR_1.0-11 broom_1.0.6 igraph_2.0.3
## [136] bslib_0.7.0